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Modern earthwork compaction rollers collect location and compaction infor- 
mation as they traverse a compaction site. These data are indirectly observed 
through non-linear measurement operators, inherently multivariate with com- 
plex correlation structures, and collected in huge quantities. The nature of 
such data was investigated at a large, atypically compacted test bed in Florida, 
USA. Exploratory analysis of this data through detrending and empirical semi- 
variogram estimation is performed. A second analysis using a sequential, spa- 
tial backfitting algorithm is used to investigate the importance of driving di- 
rection of the roller. 
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1 Modern Earthwork Compaction 

Modern compaction rollers monitor soil properties by observing stiffness characteristics 

of the soil. A vibrating drum traverses the compaction site at approximately Im/s, 

compacting approximately 20cm of material at a time. Common construction practice is 

to compact several layers of material during the construction of a new road. Each layer 

is compacted in several passes of the roller until sufficient compaction is achieved. 
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Typical construction practice is to compact in segments of road 10-15m wide and 
50-lOOm long. The roller traverses the compaction site in a snaking motion of several 



adjacent lanes. In practice, there is very little overlap between lanes (Mooney et al. 



2010). See Figure [T] for a typical compaction roller manufactured by Ammann. 




Figure 1: Ammann roller at work. 



1.1 Roller Measurement Values (RMVs) 

A typical smooth drum has a diameter of approximately Im and is approximately 2m long. 
An on-board sensor and GPS system record measurements that are together termed the 
roller measurement value (RMV). An individual RMV is an aggregate measure of a bulb 



of soil extending to a depth of approximately Im with a diameter of 0.5-0. 6m (Facas 



2009j). 

The physical nature of driving the roller down a lane with its vibrating drum causes 
other vibrational "wobbling" that remains fairly uniform over the course of the entire 
lane. Any bias this action produces will therefore be uniform over the entire lane. When 
the roller turns around and makes another pass down a different lane, the "wobbling" 
effect may be different though. This will lead to a change in the bias in the transverse 
direction, but the driving direction should remain unchanged as the new bias will be 
uniform over that entire lane. This is a cause of potential measurement error found only 
in the transverse direction. 



1.2 Florida Test Bed Data 

For a detailed investigation of roller properties, statistical characteristics, etc., a test bed 
with atypical dimensions was atypically densely compacted. A compaction roller traversed 
the compaction site in both the x- and y-directions. This Florida dataset consists of 19,145 
observations of x- and y-coordinates, soil stiffness [kg), and lane number in the x-direction 
driving and 19,975 observations in the ?/-direction driving. This analysis focuses on the 
driving direction. 

The roller first traversed the compaction site in the x-direction in a snaking fash- 
ion, first left-to-right and then back again right-to- left. The roller then traversed the 
compaction site a second time in the y-direction. The physical limitations of the site 
prohibited a snaking traversal in the y-direction, so the roller moved from bottom to top 
only. There are 29 lanes in the x-direction and 27 lanes in the y-direction. Figure |2] is 
a plot of the RMVs in the x-driving direction and the y-driving direction. Blue values 
represent high stiffness and red values represent low. An optimally compacted site would 
be uniformly blue. 
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Figure 2: Data from the test bed in Florida, USA. RMVs collected from driving in the 
x-direction (left) and from driving in the y-direction (right) are depicted. 
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2 Exploratory Data Analysis 



For this analysis, the x- and y-direction driving data are treated as two separate datasets. 
First, empirical semivariograms of the raw data were calculated using a subsample for 
computational reasons. These semivariograms exhibit aspects of non-stationarity. See 
Figure |3] for representative empirical semivariograms of both driving directions. 
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Figure 3: Empirical semivariograms of raw RMV data for x-direction driving (black) and 
?/-direction driving (red). 



2.1 Detrending the Data 

The raw data exhibits a mean trend that must be removed as a constant mean is required 
to attain second-order stationarity. By detrending the data, we can remove the mean 
trend and proceed with the analysis utilizing a second-order stationary spatial process as 
a model. 
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2.1.1 Small and Large Scale Variation 

Often times spatial data is modeled as 



y{s) = fi{s) + ais)+6{s), (1) 

where bmu{s) is the mean structure of the process, a{s) is the stochastic dependence 
structure of the process, and 6{s) is the measurement error. The mean structure is termed 
large scale variability and the dependence structure is termed small scale variability. What 
is termed mean structure and what is termed covariance structure is largely discretionary 



(Cressie 1993). 



2.1.2 Detrending Methods 

Assuming model Q, we desire a second-order stationary process a{s). Therefore, the 
data detrending process should leave some structure in the data or all that will be left is 
the noise process s{s), which is assumed uncorrelated. We have an assumption of spatial 
correlation. 

The natural first choice for detrending is fitting a linear model: /u(s) = X/3. The 
residuals of the linear model fit can then be used to estimate the semivariogram of the 
stochastic structure terms a{s) + e{s). The detrending process used included all cross 
products of the x- and y- coordinates. That is, for a 4th order polynomial, all products 
of X and y with a combined power of 4 or less were used. Empirical semivariograms were 
then calculated on the residuals of the linear model. Using a polynomial detrending of 
a 5th power generates empirical semivariograms with qualitatively identifiable nugget, 
partial sill and range parameters. This degree of detrending is desirable as all spatial 
variation is not lost and a constant mean of the residuals has been attained. 

A practical, physical explanation of the linear model parameters is not of importance. 
The goal of detrending is establishing a constant mean of the residuals, and interpretability 
of the model parameters is insignificant. 
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An alternative to linear regression for data detrending is to detrend the data using 
a nonparametric function. For this analysis, the implementation of local polynomial 



smoothing known as locally weighted scatterplot smoothing (loess) was used (Cleveland 



1979). The loess smoothing approach is based on a moving window. A polynomial is fit 
to the data in a window using robust methods. The fitted value is then the predicted 
response at the middle of the window. The window is then slid over the range of the data, 



repeating the fitting process as the window moves (Faraway, 2006). 

For this analysis, a span of 0.5 was used to reproduce empirical semivariograms like 
those of the polynomial detrend. This span corresponds to an estimated number of 
parameters of 13.5. This is approximately equivalent to a polynomial fit of 4th order, 
making this method comparable to that of a polynomial detrending. 

2.2 Fitting to a Model 

The empirical semivariograms calculated from the loess detrended data were then fitted 
to a spherical model with Cressie weights using the variof it function in R. The spherical 
model was chosen as the empirical semivariograms seemed to exhibit a linear behavior 
near the origin. The spherical model also induces sparse matrix structures, helpful for 
computation. The spherical model is defined as 



.^ , ^o(l - 1.5(V^^i) + .5(VW for he[Q,e{) 
C[n;0) = \ , [2) 

for h> 01 



where Oq is the (partial) sill and 6i is the range of the spatial process. 

Cressie weights were chosen because they are the most commonly used weights for 
fitting empirical semivariograms to a covariance model. Weighted least squares and gener- 
alized least squares require knowing the covariance structure of the semivariogram. While 



this is possible, it is hard to implement. Cressie (1985) proposed a weighting structure 



that is a compromise of weighted least squares that is no more difficult to compute than 
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ordinary least squares. 



2.3 Semivariogram Uncertainty 

Estimates of total sill, range, and nugget have very large confidence intervals. As the 



lag distance increases, the confidence interval for the total sill also increases (Nordman 



and Caragea, 2008). A simple simulation of several random fields with semivariogram 



parameters chosen to match those of the empirical semivariograms from this study was 
performed. From these random fields, empirical semivariograms were then calculated. A 
mean and standard deviation of these semivariograms was then calculated and these were 
used to calculate pointwise confidence intervals. The estimated confidence bound of the 
semivariogram starts very small for a lag distance of zero and begins expanding for larger 
lag distances. This expansion continues for larger lag distances. Decreasing the confidence 
to 75% does very little to improve the width of the estimated confidence bounds for large 
lag distances, see Figure |4j 
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Figure 4: Mean of simulated empirical semivariograms (solid line) and 95% (dashed line) 
and 75% confidence bands (dotted line). The true spherical semivariogram is depicted in 
green. 



7 



o 



o 



o 

CO 



o 

OJ 




o - 



o - 







5 



10 



15 



20 







5 



10 



15 



20 



lag dist 



lag dIst 



Figure 5: Directional empirical semivariograms and fitted spherical models from three 
polynomial detrended subsamples of RMVs of the x-driving direction subset (left) and of 
the ^/-driving direction subset (right). Dashed lines indicated fitted models, x-directional 
semivariograms are in black and green and ?/-directional semivariograms are in red and 
blue. 

2.4 Sampling Concerns 

To maintain computational efficiency, the data was subsampled for empirical semivari- 
ogram estimation. 10,000 data points were sampled from each of x- and y-direction driving 
datasets. A loess detrending of each sample was performed. This produced two detrended 
datasets from which subsamples of 2500, 3500, and 4500 data points were drawn. Direc- 
tional empirical semivariograms were then calculated in the x- and ?/-direction to generate 
a total of twelve empirical semivariograms. These empirical semivariograms were then fit 
to a spherical model. 

There was no discernible difference between the empirical semivariograms within each 
dataset. Figure [5] depicts the empirical and fitted directional semivariograms of the x- 
driving direction dataset (left) and y-driving direction (right). Since the sampled direc- 
tional empirical semivariograms are essentially identical within each dataset, we concluded 
the subsampling was adequate, i.e. the subsamphng produced a representative sample 
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2.5 Results 

The next step is a qualitative analysis of the characteristic semivariogram features. Semi- 
variograms for both driving directions exhibit similar features. For both driving directions, 
the y-directional semivariograms have a range of 0-5 and the x-directional semivariograms 
exhibit a range of 9-15. 

For the x-driving direction, the total sill for y-directional semivariograms is 22-25, 
and 15-18 for x-directional. In the y-driving direction, the total sill is 33-36 for y- 
directional and 26-30 for x-directional semivariograms. Similarly, for x-driving direction, 
the nugget for y-directional semivariograms is 5-8 and 2-5 for x-directional. For the 
y-driving direction, the nugget for y-directional semivariograms is 10-15 and 10-12 for 
x-directional semivariograms, see Table [T]and Figure [5] 

Table 1: Directional semivariogram parameters 





x-driving 


y-driving 




x-directional y-directional 


x-directional y-directional 


range 


9-15 0-5 


9-15 0-5 


total sill 


15-18 22-25 


26-30 33-36 


nugget 


2-5 5-8 


10-12 10-15 



3 Anisotropy Concerns 

Based on these observations, it is fairly safe to assume that there is no sill or nugget 
anisotropy. There does appear to be a range anisotropy between the x-directional semi- 
variograms and the y-directional semivariograms. The ratio of the range in the x-direction 
vs. the ?/-direction is approximately 5:1. The empirical semivariograms indicate a geo- 
metric range anisotropy that can be dealt with by a simple transformation of the data 



locations, (Zimmerman 1993) 



This geometric range anisotropy can possibly be explained by the compaction pro- 
cess. As the roller traverses the compaction site, it collects data every 10cm in the 
driving direction. Data is collection in the direction perpendicular to the driving direc- 
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tion approximately every l-2m. The vibrating drum is approximately 2m wide, thus the 
y-directional location of observations in adjacent lanes is l-2m apart, dependent on the 
placement of the GPS unit. Also, material is brought into the compaction site via a dump 
truck and laid down in sections. It is unclear if discontinuities exist on the boundaries of 
these sections. If they do exist, they could contribute to range anisotropy. 

Due to the nature of the driving process, data points are much more closely spaced in 
the driving direction than they are in the transverse direction. This leads to difficulties 
estimating the nugget in the transverse direction as the smallest lag distance is on the 
order of l-2m. The nugget anisotropy could therefore be explained by a vertical shift of 
the entire semivariogram caused by a measurement error in the transverse direction. This 
would essentially be a nugget anisotropy model. 

Let the true compaction process be denoted by Z{x,y) and the data we collect be 
denoted by Y{x,y) = Z{x,y) + e{y), where e{y) is a measurement error seen only in the 
y-direction. Then, the semivariogram of the Y process is jyihxjhy) = Va.T{Y{x,y) — 
Y{x + hx,y + hy)) = Var(Z(x, y) - Z{x + hx,y + hy) + e{y) - e{y + hy)) = -fzih^, hy) + 
7e(/iy). The x-directional semivariogram is then •jxih^) = 7y(/ia;,0) = 7z(^x,0) and the 
y-directional semivariogram is lyihy) = 7^(0, hy)+'^^{hy). Thus the transverse directional 
semivariogram is shifted up by the measurement error e. 



4 Driving Direction Investigation 

We utilize a state-space formulation to handle unique observation locations. Assume the 
RMVs can be decomposed into an underlying mean trend dependent on spatial location, 
driving direction, speed, and vibration amplitude, and a Gaussian spatial random process, 
(i.e. w = 1^(3 + ck), where the domain of ?/; is a lattice. Here, X is a full rank matrix 
of the fixed effects covariates and ex represents an unknown, spatially varying random 
process. The observed locations of the RMVs are then mapped to the lattice. 



Implementing a sequential, spatial mixed-effects model Heersink and Furrer (2013) 
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we can model the Florida dataset as: 



Zx = Hx^xPx + i^x0^x + ^x 

Zy Hy^i.yf3y CHyOLx ~\~ HyCKj^ ~|~ £y^ 



where cxx and cxy correspond to random variation of the layer of material being compacted 
during driving in the x- and y-direction, Hj. and H^, are the operators mapping lattice 
points to observed locations, and Sx and Sy represent the measurement error of the sensor. 
For this analysis, the lattice chosen is an 80 x 80 grid of points equally spaced on [25, 55] x 
[—.5, 33]. The size of the grid was chosen to encompass all observation locations. 

We also utilize a range anisotropy parameter p, given the empirical semivariograms 



calculated in Section 2.4 The range anisotropy is handled with a transformation of the 
coordinates. Thus p is the ratio of the range in the x-direction to that in the y-direction 
and the transformation matrix A is defined as A = diag(l,p). 



As detailed in Heersink and Furrer (2013), any additive term that can be estimated 



in a mathematically equivalent way as universal kriging can also be included in such a 
model. Splines are such an additive component that has this mathematical equivalency. 
The literature on splines is extensive and computational feasibility can be maintained, i.e. 



Wahba (1990), Filers et al. (1996), Marx and Filers (1998), Filers and Marx (2004). 



Since there was not a new layer of material added to the compaction site after com- 
pacting in the x-direction, the measurements in the y-direction are measurements of the 
same process as those in the x-direction. Thus, we should expect c — )■ 1 and either ck^, — )■ 
or OLy —7- 7y, where represents a spatially varying process only in the y-direction, e.g. a 
process representing the nugget anisotropy discussed in Section |3j Thus we would expect 
to see an empirical semivariogram of cXy to either have a sill of zero or a very small range. 
Due to measurement errors, a pure nugget model is not expected. 



The Sequential Backfitting Algorithm from Heersink and Furrer (2013) was applied 



to the data, setting c = 1, with p = 2, X is the fixed effects matrix containing all 5th 
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Figure 6: Fitted x-directional semivariograms for cx^ + Sx (black) and cx.y + £y (green) and 
y-directional semivariograms for for cx^ + Sx (red) and cxy + £y (blue) for c = 1 (left) and 
c = (right). The c = plot reproduces the curves from Figure Isl as would be expected. 



degree and lower polynomial combinations of the centered and scaled x- and y-direction 
coordinates of the roller. To create sparse matrix structures and aid in computation a 
spherical covariance function was assumed, see equation ([2]). 

The semivariogram estimation done in this study is directional. The empirical semivar- 
iograms were calculated in the driving direction. Thus, for cx^ empirical semivariograms 
were calculated in the x-direction and in the ?/-direction for cxy. 



5 Backfitting Results 

The estimated covariance parameters for otx are Ox = (12.99, 7.72)""" and 6y = (25.08, 0.39)""" 
for a.y, the estimated variances of £x and £y are = 2.68 and = 0.75, see Figure [6j 
The backfitting procedure thus reproduces the empirical x-directional semivariogram from 
the standard detrending approach. 

The range of the ay process is relatively small, thus there is no evidence to reject the 
assumption that cXy — )■ 7 from this analysis. This backfitting analysis thus reconfirms the 
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existence of a nugget effect in the ?/-direction. This would imply the "static" rolling done 
by the roller after compaction was completed is generally truly static and the material is 
not being actively compacted during this phase of construction. 

The backfitting procedure was also run for c = 0. As can be seen in the right plot 
of Figure |6} the calculated semivariograms are reproductions of the standard detrending 
approach of Section 2.1.2 and the empirical semivariograms found in Figure |5l 
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